Parameter estimation for the discretely observed 
fractional Ornstein-Uhlenbeck process and the 

Yuima R package 

Alexandre Brouste* 
Laboratoire Manceau de Mathematiques 
Universite du Maine 
Avenue Olivier Messiaen - 72100 Le Mans, France 

Stefano M. Iacus^ 
Department of Economics, Business and Statistics 
University of Milan 
Via Conservatorio, 7 - 20122 Milan, Italy 



Abstract 

This paper proposes consistent and asymptotically Gaussian estimators 
for the parameters A, a and H of the discretely observed fractional Ornstein- 
Uhlenbeck process solution of the stochastic differential equation dYt = 
—XYtdt + adW^ 1 , where (W^,t > 0) is the fractional Brownian motion. 
For the estimation of the drift A, the results are obtained only in the case 
when \ < H < §. This paper also provides ready-to-use software for the R 
statistical environment based on the YUIMA package. 

1 Introduction 

Statistical inference for parameters of ergodic diffusion processes observed on dis- 
crete increasing grid have been much studied. Local asymptotic normality (LAN) 
property of the likelihoods have been shown in [10] for elliptic ergodic diffusion, 
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under proper conditions for the drift and the diffusion coefficient, and a mesh 
satisfying 

A N — > and NA N — > +00 

when the size of the sample N grows to infinity. Estimation procedure have been 
studied by many authors, mainly in the one- dimensional case (see, for instance, 
[HI [H] and [26] in the multidimensional setting). All estimators in the previous 
works are based on contrasts (for contrasts framework, see [9]), assuming in the 
general case, that for some p > 1, as n — > +00, NA P N — )• 0. In particular, 
for Ornstein-Uhlenbeck process, transitions densities are known, and all have been 
treated, as remarked in [15] . 

In the fractional case, we consider the fraction Ornstein-Uhlenbeck process 
(fOU), the solution of 

dY t = -XY t dt + adW t H 

where W H = (W^,t > 0) is a normalized fractional Brownian motion (fBM), i.e. 
the zero mean Gaussian processes with covariance function 

EWf W t H = l - (\s\ 2H + \t\ 2H -\t- s\ 2H ) 

with Hurst exponent H e (0,1). The fOU process is neither Markovian nor a 
semimartingale for H 7^ | but remains Gaussian and ergodic (see [5]). For H > 
|, it even presents the long-range dependance property that makes it useful for 
different applications in biology, physics, ethernet traffic or finance. 

Statistical large sample properties of Maximum Likelihood Estimator of the 
drift parameter in the continuous observations case have been treated in [U HI EJ [15] 
for different applications. Moreover, asymptotical properties of the Least Squares 
Estimator have been studied in [IT] . 

In the discrete case and fractional case, we can cite few works on the topic. On 
the one hand, very recent works give methods to estimate the drift A by contrast 
procedure pH EJ] or the drift A and the diffusion coefficient a with discretization 
procedure of integral transform [25]. In these papers, the Hurst exponent is sup- 
posed to be known and only consistency is obtained. On the other hand, methods 
to estimate the Hurst exponent H and the diffusion coefficient are presented in [3] 
with classical order 2 variations convolution filters. 

To the best of our knowledge, nothing have been done, to have a complete 
estimation procedure that could estimate all Hurst exponent, diffusion coefficient 
and drift parameter with central limit theorems and this is the gap we fill in this 
paper. Moreover, estimates of H, a and A presented in this paper slightly differ 
from all those studied previously. 

In Section [2] we review the basic facts of stochastic differential equations driven 
by the fractional Brownian motion and we introduce the basic notations and as- 
sumptions. Section [3] presents consistent and asymptotically Gaussian estimators 
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of the parameters of the fractional Ornstein-Uhlenbeck process from discrete ob- 
servations. In Section [4] we present ready-to-use software for the R statistical 
environment which allows the user to simulate and estimate the parameters of the 
fOU process. We further present Monte-Carlo experiments to test the performance 
of the estimators under different sampling conditions. 

2 Model specification 

Let X = (Y t ,t > 0) be a fractional Ornstein-Uhlenbeck process (fOU), i.e. the 
solution of t 

Y t = y -X [ Y s ds + aW t H , t > 0, Y = y Q , (1) 
Jo 

where unknown parameter •§ = (X,a,H) belongs to an open subset of (0,A) x 
[g_,a] x (0,1), < A < +oo, < a < a < +oo and W H = (W t H ,t > 0) is a 
standard fractional Brownian motion [TBI HH] of Hurst parameter H e (0, 1), i.e. 
a Gaussian centered process of covariance function 

EW t H W s H = l - {t 2H + s 2H -\t- s\ 2H ) . 

It is worth emphasizing that in the case H = |, is the classical Wiener 
process The fOU process is neither Markovian nor a semimartingale for H ^ | 
but remains Gaussian and ergodic. For H > |, it even presents the long-range 
dependance property (see (5j). 

The present work exposes an estimation procedure for estimating all three 
components of $ given the regular discretization of the sample path Y T = (Y t ,0 < 
t < T), precisely 

(X n :=F nAjv ,n = 0...iV) , 

where T = Tn = NAn — > +oo and — > as N — > +oo. 

In the following, convergences — >, — and stand respectively for the con- 
vergence in law, the convergence in probability and the almost-sure convergence. 

3 Estimation procedure 

Contrary to the previous works on the subject, we consider here the problem 
of estimation of H, o and A when all parameters are unknown, using discrete 
observations from the fractional Ornstein-Uhlenbeck process. Due to the fact that 
one can estimate H and a without the knowledge of A, our approach consists 
naturally in a two step procedure. 
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3.1 Estimation of the Hurst exponent H and the diffusion 
coefficient a with quadratic generalized variations 

The key point of this paper is that the Hurst exponent H and the diffusion coef- 
ficient a can be estimated without estimating A. 

Let a = (a , . . . , ax) be a discrete filter of length K + 1, K 6 N, and of order 
L > 1, K > L, i.e. 

K K 

a k k E = for < t < L - 1 and ^ a k k L ^ 0. (2) 

k=0 k=0 

Let it be normalized with 

K 

J>l)^a fc = l. (3) 

k=0 

In the following, we will also consider dilatated filter a 2 associated to a defined by 

ak ' [ik = 2k ' for 0<k<2K. 
(J smon. 

2K K 

Since ^a 2 k k r = 2 r ^2k r a k , filter a 2 as the same order than a. We denote by 

k=0 k=0 

N-K / K x 2 



vn,* = yi yi a kXi- 



^i+k 

i=Q \fc=0 

the generalized quadratic variations associated to the filter a (see for instance 
and, finally, 

H N = - log 2 — — 

and 

i 

2 



N.a 



Theorem 1. Let & be a filter of order L > 2. Then, both estimators Hn and <7jv 
are strongly consistent, i.e. 

(Hnj&n) — > (H,a) as N — > +oo. 
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Moreover, we have asymptotical normality property, i.e. as N — >• +00, for all 

He (0,1), 



N(H N -H) AjV(0,ri(#,a)) 



and 



iV /• 

" a) AjV(0,r 2 (tf,a)) 



logiV 

where T^t?, a) and r 2 ('#, a) symmetric definite positive matrices depending on a, 
H, A and tae /liter a. 

Proof. The solution of Q can be explicited 

F t = x e- xt + a [ e- x{t ~ s) dW s H , 
Jo 

where the integral is defined as a Riemann-Stieljes pathwise integral. Let us con- 
sider the stationary centered Gaussian solution 

Y} = I e-^-^dW^. 
J —00 

We have also, 

Y} -Y t = e~ Xt (yj - y ) ^ 0. 
It is known (see Lemma 2.1]) that 

EYj(Yj-Y*) = -a 2 H(2H - l)e~ xt J° e Xu (^j\ Xv {v - u) 2H - 2 dv^j du. 

4 



Let v(t) denote the variogram of Y t \ We now show that 

Y 



«(t) = E (Y t) 2 - EF/r o t = ^\t\ 2H + r(t) 
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where r(t) = o(\t\ 2H ) as t tends to zero. Indeed, 



v(t) 



Thus, 

dv 



with 



-a 2 H(2H - 1) J° e Xu QT e~ x W(v - u) 2H ~ 2 dv^j du 
-a 2 H{2H-l) y° e Xu (^j\- Xr {t-r-u) 2H - 2 dr^jdu 



e- x{r+u \t-r + u) 2H - 2 drdu 

nu-\-t 
e~ Xw {t-w + 2u) 2H ~ 2 dwdu 

POO / PW \ 

-a 2 H(2H -I) I e~ Xw I (t - w + 2u) 2H - 2 du dw 

J0 \J max(0,w-t) J 

1 poo / pt+w \ 

—o 2 H(2H - 1) / e~ Xw ( / x 2H ~ 2 dx dw. 

2 Jo \J\t-w\ J 



i<7 2 /f(2/f - 1) r e - xw ((* + to) 2 " -2 - it - »i 2 "- 2 ) ^ 

-<7 2 F(2# - l)t 2H - 1 r e~ xt y ((1 + y) 2H - 2 - |1 - y| 2 ^ 2 ) 
2 Jo 

^ 2 #(2# - l)^- 1 ((1 + y) 2 ^ 2 - |1 - yl 2 ^ 2 ) dy +r(t) 



2 
1 
2 

cx 2 Ht 2H - l +r(t) 



_ 

<oo 

2 TJJ.2H-1 



(t) = -l* 2 H(2H - I); 2 "" 1 f; f° ((1 + y) 2 ^ 2 - |1 - yl 2 ^ 2 ) dy. 



i=i 

Therefore, we proved that 



2 

f(*) = y|t| 2H + r(*)- 



Now, applying results in Theorem 3(i)], the proof of Theorem[T]is complete 
because the following conditions are fulfilled: 

• firstly, r(t) = o(\t\ 2H ) as t tends to zero, 



• secondly, for classical generalized quadratic variations or order L > 2 (for 
instance L = 2), 

|r (4) (t)| < G\t\ 2H+1 - e - 4 
with2#+l-£ > 2#and4 > 2H+1-E+1/2 for e < 1 and any H G (1/2,1). 

□ 

Remark 1. We have two useful examples of filters. Classical filters of order L > 1 
are defined by 

f K \ (-\) l - k K\ 
ak = CL > k = ^2^{k)= 2* k\(K-k)\ P ° Ur °^ k ^ K - 

Daubechies filters of even order can also be considered (see J7j/j, for instance 
the order 2 Daubechies ' filter: 

-^=(•4829629131445341, -.8365163037378077, .2241438680420134, .1294095225512603). 
v2 

Remark 2. For classical order 1 quadratic variations (L = 1 ) and a = ( — |, 1) 
we can also obtain consistency for any value of H , but the central limit theorem 
holds only forH<\ (see 

3.2 Estimation of the drift parameter A when both H and 
a are unknown 

From |llj . we know the following result 

lim var(Ft) = lim \ [ Y t 2 dt = - 1 - -■ fi 2 ■ 
This gives a natural plug-in estimator of A, namely 



A 



2/2 2 



N 



a 2 N T (2H N + 1 
where //2,at is the empirical moment of order 2, i.e 

1 N 

-Yx 



n=l 
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Theorem 2. Let H e (|, |j and a mes/i satisfying the condition NA P N — >■ 0, 
p > 1, as N — >■ +oo. Then, as N — > +oo ; 

and 

v^fXv-A) AjV(o,r 3 (i?)), 

w/iere r 3 (#) = A (§§) 2 and 

Proof. Let us note T/v = NA^. It had been shown in [TT] that, as T/v — >■ +oo (or 

as N — > +oo), 

X 2 dt ^ k h \- 2H (5) 



^iV Jo 

and, with straightforward calculus, 



^(^rf^ X 2 dt - K H \- 2H ^j ^ JV(0, (a H n H ) 2 X- 411 - 1 ) (6) 

where kh = cr 2 r ^ +1 ^ and an is defined by Q. Let us denote /^.iv the discretiza- 
tion of the integral 



71=1 

Then 



1 r 

j"2,iV = and ^2 = «HA 2 " 



1 ' —2 7. \ , . /^r~ / 1 / ^2, 



VT N (j2 2 ,n - Ha) = V t n \ fi2,N - 7^r J X t dtj+^/T N y— J X t dt- /i 2 

As (Xt,t > 0) is a Gaussian process and Holder of order | < H < |, we have 
V^n ($2,n ~ ^ Io N X ? dt ) ^ OasiV — > +oo provided that NA P N — ► 0, 
p > 1, (see [TU Lemma 8]), we deduce from ([5]) and (J6j) that 

VT^ (h,N (0, (<x^) 2 A- 4 ^ 1 ) . (7) 

Let us introduce the following two quanitites 











M N = | 


H N 


and m = 


I?) 











Finally, results obtained in Theorem [T] and the convergence in (JTJ) gives consistency 

et ui 

2/x 2 



p 

of Mat, i.e. — > m as as JV — > +00. Let us further define 



a 2 T(2H + l] 



The derivatives of g with respect to a, H and \i2 are bounded when < A < +00 
< a < a < +00 and | < H < |. Therefore, as A^r (logA^) 2 — )• OasiV — ► +00 
we can obtain by Taylor expansion that 



'T N (g(M N ) - g{m)) A AT(0,^ 2 (m) 2 ( W *) 2 A" 4 "- 1 ' 



or 



^(aat-a) AjV(o,r 3 W) 



where r 3 (^) = g^ 2 (m) 2 (a H K H ) 2 X^ 1 = A (||) 2 , ^(.) is the derivative of g 
with respect to \i% and 

A^^A 



as — > +00. 



□ 



Remark 3. The different conditions on An raise the question of whether such a 
rate actually exists. One possible mesh is An = ^k—- 

Remark 4. As in the classical case H = \, the limit variance ^(i?) does not 
depend on the diffusion coefficient o. Let us also notice that the quantity a 2 H 
appearing in r 3 (-#) is an increasing function of H . 



4 Statistical software and Monte- Carlo analysis 

In this section we present a brief introduction to the yuima package for R statistical 
environment [21J. The yuima package is a comprehensive framework, based on the 
S4 system of classes and methods, which allows for the description of solutions of 
stochastic differential equations. Although we cannot give details here, the user 
can specify a stochastic differential equation of the form 

dX t = b(t, X t )dt + a(t, X t )dW t H + c(t, X t )Z t 

where the coefficients &(•>•), cr(-, •) and c(-, •) are entirely specified by the user, 
even in parametric form; (Z t ,t > 0) is a Levy process (for more information on 
Levy processes, see [2], [22] and (W t ,t > 0) is a fractional Brownian motion (recall 
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1 

that (W t 2 ,t > 0) is the standard Brownian motion). The Levy process (Z t ,t > 0) 
and the fractional Brownian motion (W t H ,t > 0) can be present at the same time 
only when H = |, but all other combinations are possible. The yuima package 
provides the user, not only the simulation part, but also several parametric and 
non-parametric estimation procedures. In the next section we present an example 
of use only for simulation and estimation of the fractional Ornstein-Uhlenbeck 
process considered in this paper. 

To test the performance of the estimators for finite samples, we run a Monte- 
Carlo analysis. We consider different setup for the parameters even outside the 
region | < H < | and different sample size with large and small values of T in 
order to test the performance of the estimator of the drift parameter when the 
stationarity is not reached by the process. All numerical experiments presented in 
the following have been done with the yuima package [23] . 

4.1 Example of numerical simulation and estimation of the 
fOU process with the yuima package 

With the yuima package the fractional Gaussian noise is simulated with the Wood 
and Chan method [21] or other techniques. We present below how to simulate one 
sample path of the fractional Ornstein-Uhlenbeck process with Euler-Maruyama 
method. For instance, loading the package with 

library (yuima) 

we can simulate a regularly sampled path of the following model 

X t = l-2 [ X t dt + dW t H , H = 0.7, 
Jo 

with 

samp <-setSampling(Terminal=100, n=10000) 

mod <- setModel(drift="-2*x" , dif f usion=" 1" ,hurst=0 . 7) 

ou <- setYuima(model=mod, sampling=samp) 

fou <- simulate (ou, xinit=l) 

The estimation procedure of the Hurst parameter have been implemented in qgv 
function. In order to estimate only the parameter H, one can use 

qgv (fou) 

that works also for non linear fractional diffusions (see [IS])- The procedure for 
joint estimation of the Hurst exponent H, diffusion coefficient a and drift param- 
eter A is called lse(,frac=TRUE). So for example, in order to estimate the three 
different parameters H, A and a, one can use 
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lse(fou,frac=TRUE) 



which uses by default the order 2 Daubechies filter (see Remark [T]) if the user does 
not specify the filter argument. 



4.2 Performance of the Hurst parameter and diffusion co- 
efficient estimation 

In this first simulation part, we present mean average value s and standard deviation 

for the definitions) with 



3.1 



values for both estimators #v and &n (see Section 
500 Monte-Carlo replications. This have been done for different Hurst exponents 
H and different diffusion coefficients a in the model Q, the parameter A being 
fixed equal to 2. The results are presented in Table [T] and Table [2] for different 
values of the horizon time T/v and the sample size N. 



H 


if = 0.5 


# = 0.7 


H = 0.9 




(7 


H = 0.5 


H = 0.7 


H = 0.9 


a = 1 


0.499 


0.697 


0.898 




a = 1 


1.024 


1.016 


1.081 




(0.035) 


(0.033) 


(0.031) 






(0.262) 


(0.282) 


(0.437) 


a = 2 


0.498 


0.700 


0.898 




a = 2 


2.035 


2.073 


2.213 




(0.033) 


(0.034) 


(0.033) 






(0.510) 


(0.564) 


(1.110) 



Table 1: Mean average (and standard deviation in parenthesis) of 500 Monte-Carlo 
simulations for the estimation of H (left) and a (right) for different cases. Here 
T = 100, N = 1000 and A = 2. 



H 


# = 0.5 


# = 0.7 


# = 0.9 




a 


H = 0.5 


# = 0.7 


# = 0.9 


a = 1 


0.500 


0.700 


0.900 




a = \ 


1.000 


1.001 


0.999 




(0.003) 


(0.003) 


(0.003) 






(0.025) 


(0.026) 


(0.036) 


(7 = 2 


0.500 


0.700 


0.900 




a = 2 


2.001 


2.002 


1.997 




(0.004) 


(0.003) 


(0.003) 






(0.053) 


(0.053) 


(0.073) 



Table 2: Mean average (and standard deviation in parenthesis) of 500 Monte-Carlo 
simulations for the estimation of H (left) and a (right) for different cases, and for 
T N = 100, N = 100000 and A = 2. 



Contrary to the estimation of the drift (see Section 4.3), we have consistent 
estimates of H and a for any values of T/y. Only the size of the sample N have 
influence on the performance of the estimate. 
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4.3 Plug-in for the estimation of drift parameter A 



In this second simulation part, we present mea n ave rage values and standard de- 
fer the definition) of the drift 



3.2 



viation values for the estimator A at (see Section 
with 500 Monte-Carlo replications. This have been done for different values of A 
and H in model Q, the diffusion coefficient a being fixed to 1 (see Remark [4]). 
The results are presented in Table [3] for different values of the horizon time Tjy 
and the sample size N. 





H = 0.5 


H = 0.6 


H = 0.7 






H = 0.5 


H = 0.6 


H = 0.7 


A = 0.5 


0.093 


0.214 


0.353 




A = 0.5 


0.476 


0.514 


0.605 




(0.037) 


(0.057) 


(0.069) 






(0.148) 


(0.166) 


(0.298) 


A = 1 


0.138 


0.276 


0.432 




A = 1 


0.906 


0.940 


1.005 




(0.052) 


(0.068) 


(0.078) 






(0.227) 


(0.238) 


(0.412) 



Table 3: Mean average (and standard deviation in parenthesis) of 500 Monte-Carlo 
simulation for the estimation of A for different values of H and A. Here a = 1 and 
T N = 1 and N = 100000 (left) and T N = 100 and N = 1000 (right). 



We can see in Table [3] that the values of Tjy is important for the estimation of 
the drift. Actually, the consistency of the estimates are valid for increasing values 
of T/v and decreasing values of the mesh size A^. Moreover, the bigger H, the 
harder the estimation of the drift parameter. This phenomena can be explained 
by the long-range dependence property of the fOU process. It is the same for A ; 
as A increases, its estimation is harder (see Remark Q. It can be explained by the 
fact that when A is bigger, the fOU process enters faster in its stationary behavior 
where it is more difficult to detect the trend. 

Finally, in order to illustrate the asymptotical normality for the estimator A of 
A, we present in Figure [T] the kernel estimation of the density. 



Acknowledgments 

We would like to thank Marina Kleptsyna for the discussions and her interest 
for this work. Computing resources have been financed by Mostapad project in 
CNRS FR 2962. This work has been supported by the project PRIN 2009 JW2STY, 
Ministero dell'Istruzione dell'Universita e della Ricerca. 



References 

[1] B. Bercu, L. Coutin, and N. Savy. Sharp large deviations for the fractional 
Ornstein-Uhlenbeck process. Teoriya Veroyatnostei i ee Primeneniya, 2010. 



12 




Figure 1: Kernel estimation for the density of ( y/T^ ( — A ] ] , M = 

\ \ / / m=l...M 

5000, for T/v = 1000 and T/v = 100000 (fill line) and the theoretical Gaussian 
density W(0, r 3 (0)) (dashed line) for = (A, a, H) = (0.3, 1, 0.7) (for the value of 
r 3 (i?) see Theorem |2|. 



[2] J. Bertoin. Levy Processes. Cambridge University Press, Cambridge, 1998. 

[3] C. Berzin and J. Leon. Estimation in models driven by fractional brownian 
motion. Annates de I'Institut Henri Poincare, 44(2):191-213, 2008. 

[4] A. Brouste and M. Kleptsyna. Asymptotic properties of MLE for partially 
observed fractional diffusion system. Statistical Inference for Stochastic Pro- 
cesses, 13(1):1-13, 2010. 

[5] P. Cheridito, H. Kawaguchi, and M. Maejima. Fractional Ornstein-Uhlenbeck 
processes. Electronic Journal of Probability, 8(3):1-14, 2003. 

[6] I. Cialenco, S. Lototsky, and J. Pospisil. Asymptotic properties of the max- 
imum likelihood estimator for stochastic parabolic equations with additive 
fractional Brownian motion. Stochastics and Dynamics, 9(2): 169-185, 2009. 

[7] I. Daubechies. Ten Lectures on Wavelets. SIAM, 1992. 

[8] D. Florens-Zmirou. Approximate discrete time schemes for statistics of diffu- 
sion processes. Statistics, 20:263-284, 1989. 

[9] V. Genon-Catalot. Maximum constrast estimation for diffusion processes from 
discrete observation. Statistics, 21:99-116, 1990. 



13 



[10] E. Gobet. Lan property for ergodic diffusions with discrete observations. 
Annales de VInstitut Henri Poincare, 38(5):711-737, 2002. 

[11] Y. Hu and D. Nualart. Parameter estimation for fractional ornstein-uhlenbeck 
processes. Statistics and Probability Letters, 80(11-12):1030-1038, 2010. 

[12] J. Istas and G. Lang. Quadratic variations and estimation of the local hlder 
index of a gaussian process. Annales de VInstitut Henri Poincare, 23(4) :407- 
436, 1997. 

[13] J. Jacod. Inference for stochastic processes. Statistics, Prepublication 683, 
2001. 

[14] M. Kessler. Estimation of an ergodic diffusion from discrete observations. 
Scandinavian Journal of Statistics, 24:211-229, 1997. 

[15] M. Kleptsyna and A. Le Breton. Statistical analysis of the fractional Ornstein- 
Uhlenbeck type process. Statistical Inference for Stochastic Processes, 5:229- 
241, 2002. 

[16] A. Kolmogorov. Winersche Spiralen und einige andere interessante Kurven in 
Hilbertschen Raum. Acad. Sci. USSR, 26:115-118, 1940. 

[17] C. Ludena. Minimum contrast estimation for fractional diffusion. Scandina- 
vian Journal of Statistics, 31:613-628, 2004. 

[18] B. Mandelbrot and J. Van Ness. Fractional Brownian motions, fractional 
noises and application. SIAM Review, 10:422-437, 1968. 

[19] D. Melichov. On estimation of the Hurst index of solutions of stochastic 
equations. PhD thesis, Vilnius Gediminas Technical University, 2011. 

[20] A. Neuenkirch and S. Tindel. A Least Square-type procedure for parameter 
estimation in stochastic differential equations with additive fractional noise. 
preprint, 2011. 

[21] R Development Core Team. R: A Language and Environment for Statistical 
Computing. R Foundation for Statistical Computing, Vienna, Austria, 2010. 
ISBN 3-900051-07-0. 

[22] K. Sato. Levy Processes and Infinitely Divisible Distributions. Cambridge 
University Press, Cambridge, 1999. 

[23] YUIMA Project Team, yuima: The YUIMA Project package (unstable ver- 
sion), 2011. R package version 0.1.1936. 



14 



[24] A. Wood and G. Chan. Simulation of stationary Gaussian processes. Journal 
of computational and graphical statistics, 3(4):409-432, 1994. 

[25] W. Xiao, W. Zhang, and W. Xu. Parameter estimation for fractional ornstein- 
uhlenbeck processes at discrete observation. Applied Mathematical Modelling, 
35:4196-4207, 2011. 

[26] N. Yoshida. Estimation for diffusion processes from discrete observations. 
Journal of Multivariate Analysis, 41:220-242, 1992. 



15 



